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Abstract 

The kinetics of coherent Cu rich precipitation in Fe-Cu and Fe-Cu-Ni alloys during thermal 
ageing have been modeled by Atomic Kinetic Monte Carlo method (AKMC). The AKMC is pa- 
rameterized by existing ab-inito data to treat vacancy mediated diffusion which is depend on local 
atomic environment. A nonlinear semi-empirical time adjusting method is proposed to rescaled the 
MC time. The combining AKMC and time adjusting method give a good agreement with experi- 
ments and other simulations, including advancement factor and the Cu cluster mobility. Simulations 
of ternary alloys reveal Ni has a temporal delay effect on Cu precipitation. This effect is caused by 
the decreasing of diffusion coefficient of Cu clusters. And the reduction effect of diffusion coefficient 
weakens with larger Cu cluster size. The simulation results can be used to explain the experimental 
phenomenon that ternary Fe-Cu-Ni alloys have higher cluster number density than corresponding 
binary alloy during coarsening stage, which is related to cluster mobility. 

PACS numbers: 81.30.Mh, 66.30.-h, 64.75.Nx, 05.10.Ln 
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I. INTRODUCTION 



Irradiation induced precipitation is believed to be one major reason for the degradation 
of mechanical properties of alloys in radiation environments. In a-Fe, which is the basics 
of ferritic steels, Cu rich precipitation from supersaturated matrix is greatly accelerated by 
irradiation. This kind of Cu rich precipitation is the primary reason of embrittlement for 
reactor pressure vessel (RPV) steels at low doses compared to the so called "late blooming 
phases" of MnNi rich precipitationsi at high doses. Experiments show third-party elements 
such as Ni, Mn and Si etc. also exist in Cu rich precipitations and may have influence on 
Cu precipitation kinetics. In those elements, Ni and Mn are the richest elements in typical 
RPV steels^. Possibly, they also have the strongest influence on Cu precipitation. Miller et 
al^i^ found that a Fe-Cu-Mn model steel has the cluster density approximately an order of 
magnitude higher than that of Fe-Cu steel, and RPV steels with high nickel content may have 
retarded precipitation growth as evidenced by smaller cluster size. Meslin et air- discovered 
that advancement of precipitation is lower in the presence of Mn and Ni, suggesting they 
may delay the copper precipitation. The common summaries of experiments are higher 
precipitation number density or lower advancement is found in ternary alloys within same 
ageing time or radiation dose of binary alloy. 

To further reveal the mechanism of the influence by third-party elements, atomic level 
computer simulation solutions seem to be attractive. Atomic kinetic Monte Carlo (AKMC) 
method based on diffusion of point defects has become an effective research tool on precipi- 
tation for having detailed information on atomic configuration in full time scale and being 
convenient to separate different factors. Vincent et air- studied the effect of Mn and Ni on 
Cu precipitation during radiation flux, the simulation results show Mn containing alloy has 
slightly smaller cluster size, while Ni seems to have little influence. Bonny et air- applied 
an artificial neural network (ANN) powered AKMC to study the precipitation of a ternary 
Fe-Cu-Ni alloy, and found the peak density of clusters increased by about 29% than binary 
alloy. However, due to the time evolution model used by previous works, the kinetics is 
not represented in view of real time evolution, e.g. MC time scale is found incomparable to 
experiments^, which limited further comparison. 

In this study, we focus on the effect of Ni on Cu precipitation kinetics. To accomplish our 
aims, an AKMC approach is applied to simulate thermal ageing of Fe-Cu and Fe-Cu-Ni 
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alloys. The following content of this paper is divided into three parts. Detail of the com- 
putational methods used in the AKMC is presented in the first part. The parameterization 
is given and a nonlinear time adjusting method based on post-processing of AKMC data is 
proposed in order to reflect the kinetics correctly. In the second part, the simulation results 
are reported. Firstly, the simulation results of precipitation kinetics of the Fe-Cu binary 
system at varying temperatures are used to verify applicability of AKMC parameters and 
combined time adjusting method. Then the results of alloys with different Ni content, aged 
at same temperature, are compared. In the final part of the paper, the effect of Ni on Cu 
precipitation kinetics is discussed. 



II. METHODS 



In the AKMC simulation, the precipitation process is induced by thermal ageing. Initially, 
Cu and Ni substitutional atoms are randomly introduced into a-Fe matrix. A single vacancy 
is randomly introduced into the system to treat vacancy mediated diffusion. 



A. AKMC simulation model 



An AKMC code has been developed at Nanjing University of Science and Technology, 
with rigid on-lattice model (RLM) and Bortz- Kalos- Lcbmvitz (BKL) algorithni^"^. When 
the vacancy lies on the first nearest neighbor lattice site of one atom, the probability of 
the position exchanging between this atom and the vacancy is obtained by the Arrhenius 
equation, 

= z/xexp ( ) , (1) 



where i>x is the attempt frequency for atom species X, which is related to local vibration 
modesi^. In the present paper, Ux values for all atom species are taken as one independent 
constant of 6x10^^ s"^, which is on the same order of Debye frequency. 

The time evolution of one KMC step is given by the summation of jump frequency of 
every possible exchanging. 

In R /^x 

(2) 

where, i? is a uniform random number between and 1. 
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B. Activation energy model 



The activation energy in eq. ([T|) plays the key role in diffusion dynamics. Comprehensive 



11 



-13 



description on energy models with heuristic formulas has been reported in Ref. 
Recently Vincent et alM made a critical review of these models. It is also possible to directly 
predict the activation energy by ANN AKMO^^"— . Here, we use the final initial system 
energy model (FISE)^iiiii^, which has the same form of the Kang- Weinberg decomposition^^. 
For the situation of atom species X exchanging with a vacancy lying on the first nearest 
site, the activation energy writes as, 

Ea = ^^^^^^^ + Qx, (3) 

where Eini and Efni are the initial and final system energies, respectively. Qx is the migration 
energy of atom species X in a-Fe matrix. 

The interaction of atoms is ranged up to second nearest neighbor and under two-body 
approximation, thus the system energy at a specific state is evaluated by summation of the 
energies of pairwise bonds, as following, 

^5 = + tt4iB, - 4'-v, (4) 

i=l j=l i=l j=l 

where i represents the bonds are counted in both first (« =1) and second {i =2) nearest 
neighbor sites. Zi is the coordination number at each distance, Aj is a neighbor atom of the 
jumping atom X , Bj is a neighbor atom of the vacancy, e is bond energy. X-V bond at first 
nearest neighbor distance is double counted, thus is subtracted in the third term. 

Basic information on atom interactions is usually obtained by ab-initio calculations in 
a multi-scaled fashion. Results obtained by Vincent et al^^^^^ has been opted in current 
study. The atomic interactions in Fe-Cu-Ni-Mn-Si system were studied using Projector 
Augmented- Wave (PAW) and UltraSoft Pseudo Potential (USPP) in their works. It was 
found that the USPP results reached a better agreement with experiments^. The values of 
Qx for Fe, Cu and Ni are 0.62eV, 0.54eV and 0.68eV respectively^^. In the present work, we 
fitted pairwise bond energies using a method similar to the description in Ref. [g]. The ah- 
initio data of cohesive energy, mixing energy, binding energy and vacancy formation energy 
are expanded by following relations in RLM, 

Et = ^e^tx + ^e^tx. (5) 



(Fe) = -441,, - 34lFe + 841x + ^4lx " " (6) 

T^Hi) (T^\ _ I (j) I (i) _ («) 

^XY \^^) — ^Fe-Fe ^Fe-X '^Fe-Y ^X-Y; \' ) 

E{r (X) = 8eil^ + Qe^tv - - S^U, (8) 

To balance the ratio between bond energies of X- Y pair on the first nearest and second 
nearest distance, extra equations are needed. Vincent et al^ assumed a constant ratio 

n 

between the second nearest bonds of X-X pair and Fe-Fe pair. Similarly, in Ref. [Ill, except 
the Cu-V bonds, the energy of a second nearest X-Y bond is half of the first nearest X-Y 
bond. Here we added extra equations by assuming the interfacial energies on {100}, {110} 
and {111} planes have negligible difference, as the following equations, 

pint{110} _ pmt{100} /QN 

^XY — ^XY ' \y) 

pint{110} _ pint{lll} 

^XY — ^XY ' \^^) 

The interfacial energies can be expanded as equations, 

^int{100} _ _2 (1) _ (2) 4 (1) 2£(2) _ 2e(l) _ Q^N 

^XY ~ ^^X-X ^X-X ' ^^X-Y ' '^'^X-Y '^^Y-Y ^Y-Y^ \^^) 

■ V2, (12) 



j-,int{110} 
^XY 



^x-x ^x-x ' ^*^x-Y ^ ^^x-y ^y-y ^y-y 



pmt{lll} 
^XY 



.V^) _3p(2) +5.(1) +6.(2) _3J2) 



/v^, (13) 



Practically, our fitting is divided into two steps. In the first step, equations ([5]) ([6]) ([7]) ([9]) ffTOl) 
are used to get the interactions between atoms. Then, in the second step, equations ([7]) and 
(IH]) are used to get the interactions between atoms and vacancy, where the interactions be- 
tween atoms obtained by first step are considered to be known. In both steps, the equation 
sets are over-determined, the Moore Penrose generalized inverse matrix method was used to 
get the least square solution. 

The fitted pairwise bond energies e^-y given in Tabled It's worthwhile to notice that 
because the fitting equation sets are over-determined, not all equations are exactly equal on 
both sides at the end. Actually, only cohesive energies of Fe, Cu, Ni were exactly fitted, 
the values are -4.28eV, -3.49eV and -4.34eV respectively. Comparison of other energies is 
presented in Table [Tll The binding energies of Cu-Cu and Cu-V are obviously smaller than 
the ah-intio data. Though, mixing energy of Cu in a-Fe matrix was fitted well. We think 
it is because Cu-X pairs have strong many-body contribution that cannot be reproduced 



TABLE I. Fitting result of pairwise bonds of atomic interactions 



Bond X- Y 







i=2 


Fe-Fe 


-0.6856 


-0.5125 


Fe-Cu 


-0.5696 


-0.4429 


Fe-Ni 


-0.7099 


-0.5201 


Cu-Cu 


-0.5610 


-0.4153 


Cu-Ni 


-0.6600 


-0.4764 


Ni-Ni 


-0.6903 


-0.5263 


Fe-V 


-0.2117 


-0.1035 


Cu-V 


-0.1853 


-0.1611 


Ni-V 


-0.2601 


-0.2767 



by pairwise bonds. As is shown in Ref. |21|, cluster expansion can get better fitting. It is 
also need to mention that for Cu and Ni, the reference structure for cohesive energies and 
mixing energies in Vincent's original work are fco^^i^^. Energy difference between fee and bcc 
is 0.036eV for Cu by Domain and Becquart^^ and O.lOeV for Ni by Mishin et alM, though 
these small differences are ignored in fitting. 

As suggested by Soisson and Fu^i, it seems necessary to introduce non-configurational 
entropy into the energy model to get better solubility of Cu, which is done by adding temper- 
ature dependent term for the Fe-Cu bond energies. We also used this method. Since saddle 
state is more related to initial state, this modification is only introduced into initial system 
energies. So Eini becomes E^^- with its pairwise bonds modified to s^^'_y = £^x-y ~ -^x-y-^- 
The A values for Fe-Cu pair are obtained by fitting Cu solubility in a-Fe, corresponding 
non-configurational entropy AS'„c equals I.S/cbj for Fe-Ni and X-V pairs, AS'„c were set as 
1.0 /cb- The A factors for related pairs are listed in Table UTTl 



C. Time adjusting 

The time adjusting has special importance for AKMC thermal ageing simulations applied 
to solid solutions. For one reason the computation cost requires relatively smaller simula- 
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TABLE II. Comparison of targ et ab-intio thermodynamic properties and fitted results in RLM 
Target values are from Ref. [l9| (except the marked ones) 



target fitted target fitted 





0.55 




0.5562 4^ , (Fe) ^ 


0.065 


0.0661 




n 1 7 

-U. 1 ( 






u.uz 


U.U/Oo 


4^«k (Fe) 


0.16 




0.1075 E'^liFe) 


-0.10 


-0.0440 


4^fk (Fe) 


0.05 




0.0421 i^S^Fe) 


-0.02 


-0.0015 


(Fe) 


2.00 




1.9658 4^;^/(Fe)l 


0.16 


0.0897 


e{;' (Cu) 1 


1.05 




1.0412 4lv(Fe)l 


0.18 


0.1272 


(Ni) 


0.60 




0.5993 E'ill (Fe) 
^5!?^(Fe) 


0.03 
0.17 


0.0241 
0.1656 


^ The original value is 1.6eV, 


sij 


^nificantly higher than ab-intio result by Soisson et al.—, also higher than 


prediction by interatomic potentials'^—. So we modified it. 
^ Fitting of original value gives incorrect number density, i.e. density decrease with increasing Ni content, 


so modified intentionally to 
Modified according to Ref. 


a 
23 


level according to the ternary potential of Bonny et al^. 






TABLE III. Entropy contribution factor A for pairwise 


bonds 




Bond X- Y 






^X-Y 










i=l 




i=2 


Fe-Cu 






8.8450 xl0"6 




6.8774x10"*^ 


Fe-Ni 






6.9518x10"*^ 




5.0930 xl0-*5 


Fe-V 






7.8816x10"*^ 




3.8534x10-6 


Cu-V 






6.5202x10"*^ 




5.6686x10-6 


Ni-V 






5.9914x10"*^ 




6.3738x10-6 



tion box which usually causing significant higher vacancy concentration and under-estimated 
time evolution. Moreover, phase separation progress results in evolving equilibrium vacancy 
concentration, while number of vacancies introduced in AKMC is fixed and integral. This 
difference cause time under-estimate degree evolves nonlinearly. We used the model de- 



scribed in Ref. 
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physical time is rescaled by the ratio between vacancy concentration 



within matrix in MC model and equilibrium vacancy concentration of matrix, 



C^c (M) 



(14) 



Equilibrium vacancy concentration of matrix is determined by vacancy formation energy in 



matrix, Cy' (M) = exp 



f )exp 



(M) 



AS* is taken as 1.0 as we mentioned above. 



Vacancy concentration within matrix in MC model is given by. 



C^^ (M) 



fM 



(15) 



NXu' 

where, Xu is the concentration of matrix atoms(Fe) in the box, A'^ is the total atom num- 
ber, is the ratio between equilibrium vacancy concentration in the solid solution and 
concentration in pure matrix. 

Our concern focuses on the /y' value. There has not been a physical analysis for this 



value. In Ref. 
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and 
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it is regarded as fraction of time spent by the vacancy in Fe 
matrix, thus in their simulations the value was computed by checking if first and second 
nearest neighbors of vacancy contain Cu atom. One difficulty of this method is that because 
the system evolves nonlinearly with time, the ideal step length is nearly impossible to be 
obtained on the fly. As a consequence, the slope of (t) curve is slightly deviated at the 
beginning of simulation, and data fluctuation continuing become remarkable stronger with 
MC time increasing, especially in the coarsening stage. Here, we propose a semi-empirical 
way to get the value. According to Lomer— , for a dilute solid solution, the value for 
ideal solution state can be determined physically, as well as the state when precipitation 



completed. Soisson and Fu^ have given the expressions for Fe-Cu system. 



(0) 



1 — ZiCcu — Z2Ccu 



[1 - ZiCcn - Z2Ccu) + Xl^iCcuGXp 

i=l 



^(oo) 



1/ 



1 



Ccu 



exp 




(16) 



(17) 



1 - Ccu " \^ A;bT 
The expressions give the upper and lower bounds for the value. Since decrease 



monotonically, an S-shape sigmoidal function is expected to link the bounds. The cluster 
number density, mean cluster size and advancement factor, also evolve with time, there 
should be a mapping relationship between the evolution of /y' and those quantities on the 
basis of time. Since the mentioned physical quantities can be calculated in post-processing, 



8 



if a sigmoidal function and mapping relationship of several special points is given, a semi- 
empirical function of evolution can be obtained. 

Since it is known evolution of advancement factor follows the Johnson-Mehl-Avrami 
(JMA) law, C,{t) = 1 — exp [— (t/r)"]. should evolve no slower than the order of JMA 
law. The sigmoidal Hill equation— has been chosen to represent the evolution, the 
formula is given as below, 

log(/)-log(/o) 1 



iog(/oo)-iog(/o) i + (w^or' 

where treai is the real time, to and p are unknown parameters. 
Use a substitution g = log (trcai), we get the following, 

log (/) - log (/o) 1 



(18) 



(19) 



log (/oo) - log (/o) 1 + exp [{g - go) / dx] ' 

where go = log (to) and dx = 1/p. 

The expression on right of eq. (fT9|) is known as Boltzmann equation. The /y' (t) calcu 



lated in Ref. 
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got perfect fitting by this equation, with go = —0.82 and dx = 0.25 (using 
10 as the base of log function). This implies the eq. ( IT9|) has the correct "evolving speed". 
So the next step to get the semi-empirical solution for /y' is to fit the unknown parameters 
of ^fQ and dx. 

We have chosen the number density of Cu clusters to give a mapping relationship of 
special points between post-processing data and /y^. The reason is that this variable does 
not evolve monotonically which making distinguish of different stages easily. On the number 
density evolution, two regions have the best analysis properties making them special. One 
is the time when nucleation starts, the other is the growth stage. When nucleation starts, 
binding of vacancy to clusters should has the most significant change, since random embryos 
from thermal fluctuation become stable nucleuses, this will cause a rapid change of the slope 
of the (t) curve. So we assumed the time nucleation starts , identified by when number 
density starts increasing, is mapped to the maximum curvature position ( = 0) of /y (t) 
curve, that is ((70 + dx -In^b — , (3 + ^/O) /6) for Boltzmann equation. Similarly, at 

growth stage, new nucleus stops to form, a maximum slope is expected since contribution 
to vacancy binding from nucleation is lost. We assumed the time growth starts mapped 
to the infiection point ( = 0), that is {go, 0.5), and identified by number density first 
reaches 90% of peak density. The mapping relationship is represented in Fig. [TJ MC time 
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position \ 
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log (treal) 



FIG. 1. Real time fitting scheme by relationship between number density and /y . 

of those two positions can be directly read from number density evolution, namely for 
time when nucleation starts and for time when growth starts. According to eq. ( fT4l) . 
the MC time can be obtained by integration as a function of real time. Thus by solving 
the problem of equations ( 120|) and (1211) with proper optimization method, (70 and dx can be 
fitted. Corresponding relationship between each MC time and real time can be obtained by 
numerical integration within an acceptable margin of error. This semi-empirical method is 
based on the post-processing of the number density data from AKMC and avoids complicate 
record step length setting and data fluctuation. The obtained j}f (t) relationship is nonlinear 
as the vacancy bind energy evolution. 

/•nuleation— start 

NXmC'^ (M) / (1//^) dWi = t^^ic, (20) 

Jo 

p growth— start 

NXmC;^ (M) / (1//^) dt,eal = t^MC ' ^MC (21) 

J nuleation— start 



III. RESULTS 



A. Cluster identification 



Cu rich clusters are identified by counting number of "bonds" linked to lattice sites. Those 
bonds are ranged up to second nearest neighbors as the activation energy model defined. 

10 



Since Fe atom is absent inside clusters in our results, which consist with other simulations, 
Fe atoms are ignored during the identification. Because of the time adjusting method we 
used, a lesser identify rule has been used, to be recognized as part of a cluster, one atom 
is required to be linked with at least three bonds. This geometrically causes the smallest 
cluster recognized is a tetrahedron. 



B. Fe Cu binary system kinetics 

We applied the model above to simulate the precipitation kinetics of a Fe-1.34 at. % Cu 
alloy (about 1.4 wt. % Cu), during thermal ageing at four different temperatures, 663K, 
713K, 773K and 873K respectively. 

The advancement factor ^, defined by eq. fl22l) , represents the completeness of precip- 
itation progress, and is a basic property which can be measured by several characteristics 
techniques. 

c _ <^Cu (0) - Ccu jt) . 

^ ^ ' ~ Ccu (0) - Ccu (oo) ' ^^'^ 
in the equation, Cqu is the concentration of solute Cu atoms. Ccu(O) is the initial Cu 
content, and Ccu(oo) is the solubility of Cu in a-Fe which is estimated by the following 
equation, 

ASZ^\ ri?6u(Fe)l ^23) 



Ccn (oo) = exp y—^j exp 

where, AS*^^" is the non-configurational entropy, and E^^ (Fe) is the mixing energy of Cu. 
To obtain the advancement factor, a simulation box containing 64x64x64 bcc unit cells 
was used. The evolution of the advancement factor is shown in Fig. O The curves have 
good agreement with experiments overall the temperature range. Exceptions are, at low 
temperatures 663K and 713K, the kinetics are slower than experiment at the beginning of 
precipitation, and at high temperatures 773K and 873K, the kinetics are faster than exper- 
iments when ^ > 0.6. Even though, several simulation tests in bigger boxes (128x128x128) 
with C, up to 0.7 show large box size seems to have better agreement. And, using more strict 
identification rule, i.e. cluster no smaller than 10 atoms, the kinetics will only be slightly 
slower at the beginning. 

Advancement factor is still not enough to refiect all aspects of precipitation kinetics. A 
simulation box of 128x128x128 unit cells has been used to get the number density and the 
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FIG. 2. Evolution of the advancement factor of Fe-1.34 at. % Cu under different temperatures. 



mean cluster radius of Cu clusters. As is shown in Fig. [3]^a), the number density and mean 
cluster size evolution has similar tendency of experiments. But the absolute value of the 
calculated number density is nearly a magnitude of order higher than the experiment data. 
In our AKMC result, after 10° h clusters are expected to have larger size than experiment 
detection limit, since the coarsening has already started. However, the number density after 
10° h is still higher than experiments. So this kind of deviate is more than identification 
difference. To our knowledge, Castin et alM gave the most satisfied result by hybrid AKMC 
method with ANN trained by an embedded atomic (EAM) potential^^. According to classical 
nucleation theory, the critical nucleation size is determined by volume free-energy change 
and interfacial energy of Cu cluster. The interfacial energies evaluated by our pairwise bonds 
are Ej^S°°^=0.5035J/m2, £;J,°S^°^=0.4142J/m2 E'^^^l^^^ =0.U65 J/m^. And molecular 



statics calculation shows the EAM potential from Ref. |26| gives 4ccr^=0-3789 J/ m 
^Sc"°^=0-4113J/m2 and Ej^S"^=0.5130J/m2. The interfacial energies evaluated by both 
methods are close. However, as we mentioned before, the binding energies of Cu atoms are 
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FIG. 3. Precipitation kinetics of Fe-1.34 at. % Cu under 773K, (a) cluster number density evolution 
Np{t), (b) mean cluster radius evolution rp{t). Experimental results are from Kampmann and 



Wagner: Ref. 



34, Goodman: Ref. 



35, Mathon and Barbu: Ref. 
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significantly under-estimated by pairwise bonds, the volume free-energy change should also 
be under-estimated. Thus the critical nucleation size will be under-estimated, result in a 
higher peak number density, as well as the density in coarsening stage. 

It was first suggested by Soisson et alM'^ that the Cu precipitation in a-Fe may favor 
a coagulation mechanism caused by highly mobile Cu clusters. This theory was recently 
confirmed by a hybrid AKMC method^^. We noticed that the difference between coagulation 
and emitting-absorbing mechanism only become marked in long term coarsening stage, our 
AKMC results seems to be not strong enough to be proved "correct". It was mentioned 
that their Object KMC part is controlled by several parameters of Cu clusters, i.e. lifetime, 
diffusion coefficients and dissolution probability. Therefore we used the AKMC simulated 
mobility of VCua? clusters at 773K. Fig. H] shows the diffusion coefficients, lifetime and 
dissolution probability versus time, each data point was simulated for tens of thousands 
of times for statistics. The diffusion coefficients and lifetime have a good agreement with 
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FIG. 4. Cluster mobility of VCuat under 773K, (a) diffusion coefficients D^v and lifetime tat versus 
cluster number (b) dissolution probability versus cluster number. Squares: this work; line: 



5-spline line of data from Ref 
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Ref. |33|. The dissolution probability by our model shows a different pattern, it does not 
decrease quickly with cluster size. However, th e dissolution probability calculated by us 
is always lower than the probability in Ref. 
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It is apparent only a system with large 



dissolution probability will tend to emitting-absorbing mechanism. According to their work, 
even if dissolution is forbidden which means a dissolution probability of 0, the evolution 
still has good agreement with slight larger cluster size compared to experiments. So our 
model actually can reproduce the coagulation mechanism, rather than emitting-absorbing 
mechanism. 

It has been shown from above results, our model reproduces consistent results with ex- 
periments and other simulations over a temperature range. The advancement factor was 
reproduced well, as well as the Cu cluster mobility. Although the cluster density is over- 
estimated by our model, it is a drawback of pairwise bonds. Nevertheless, good description 
of Cu precipitation kinetics in binary Fe-Cu system has been obtained by our AKMC model 
and combined time adjusting method, the very first step needed for the study of Ni effect. 
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C. Fe Cu Ni ternary system kinetics 



We performed thermal ageing simulations at 823K on Fe-1.34 at. % Cu at. % Ni 
alloys, the x values are 0, 0.5, 1.0, 1.5 and 2.0, respectively. The simulations were conducted 
in a box containing 128x128x128 unit cells. 

By direct visual observation from snapshot, we found that Ni appears to remain random 
distribution during the simulation, though some Ni atoms embellishes on Cu clusters like 
strawberry. Based on this fact, in the ternary alloys, one more rule was added to cluster iden- 
tification, that is clusters should have only one pure Cu core no smaller than a tetrahedron. 
And, the bounds of is extended as follows. 



E (1 - ZiCx - Z2Cx) 




(25) 

The advancement factor evolution is shown in Fig. [5]^a). Starting from about 2xl0"^h, 
the precipitation kinetics become slower in four ternary alloys. When refer to Fig. [5](b) of 
number density evolution, we can see that, before number density reaches the peak number 
density the evolution of number density in ternary alloys is even slight faster. However, 
when passed the peak number density position, the evolution of number density in ternary 
alloys become obviously slower than the binary Fe-Cu alloy. The time 2xl0"^h actually 
corresponds to the time at peak density. So from these observations it has been confirmed 
that the precipitation kinetics of ternary alloys appears to be "delayed" in the coarsening 
stage. 

At the peak density point, the ternary alloys containing 1.5 at. % and 2.0 at. % Ni have 
about 10% higher number density than the binary alloy, while the other two ternary alloys 
are about 3% higher. It was found the peak number density of Cu clusters in Fe-1.13 at. 
% Cu-1.36 at. % Ni is about 29% higher than Fe-1.13 at. % Cu by ANN AKMC in Ref. 
Q. And the experimental results by Buswell et a/.— gave about 34% higher peak number 
density by means of larger area under density-size distribution. Our result reproduced a 
similar tendency. The increasing peak number density is possibly caused by Ni effect on 
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FIG. 5. Precipitation of Fe-Cu-Ni ternary alloys under 823K, (a) advancement factor evolution 
(t) , (b) cluster number density evolution Np (t) . 



nucleation. Al-Motasem et ai— revealed CumNi„ clusters have lower formation energy than 
pure Cu clusters. Seko et alM also gave a similar result by ab-initio calculation. So the 
reason for a higher peak number density of Cu clusters is that Ni atoms act as nucleation 
centers for Cu precipitates and promote nucleation of Cu clusters. This also explains the 
more rapid evolution of ternary alloys during nucleation stage. 

The mean surface area density of Ni atoms surrounding Cu cluster is shown in Fig. [6l 
The surface area density of Ni increases quickly with time and then saturates. The more 
Ni contains in a ternary alloy, the higher saturated surface area density of Ni is formed. 
Actually, a linear relationship has been found between the saturated surface area density of 
Ni and the Ni content of the alloy, with a slope around 1.46 nm"^/(l%Ni). The figure also 
confirmed visual observation that few Ni atoms embellish on Cu clusters, since the density 
is quite low. 
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FIG. 6. The evolution of surface area density of Ni atoms on clusters (Ni) [t]. 
IV. DISCUSSION 

One obvious effect of Ni on Cu precipitation is that higher cluster number density found 
in ternary alloys during coarsening stage compared to binary Fe-Cu alloy. In other alloy 
systems, similar phenomenon caused by third-party element atoms can also be found. For 
example, in Al-Sc-Zr alloys, addition of Zr is known to have an effect of producing higher 
density of LI2 structure precipitation Al3Sc3;Zri_j; during the coarsening stage^"-^. It will 
be interesting to compare the precipitation kinetics of these two alloys system. We plot 
the mean size evolution versus advancement factor of ternary Fe-Cu-Ni alloys in Fig. El^a). 



44 and 



45 



within the 



As for Al-Sc-Zr system, the kinetic data are extracted from Ref. 
time rage from 0.3s to 0.5s, which seems comparable to our simulation time range. The 
cube root of mean size of precipitated AlsSc^^Zri.^; clusters versus total number of atoms in 
precipitations is represented in Fig. WO^)- The difference is quite remarkable, little change 
occurred by Ni addition. However, Zr addition significantly refines Al3Sci;Zri_a: precipitation. 
Clouet et al^ revealed that the diffusivity difference of Zr and Sc make Zr-rich external shell 
forming around the precipitation, which blocks Ostwald ripening. Based on the comparison 
from the figures, Ni plays a different role. The higher number density in Fe-Cu-Ni ternary 
alloys during coarsening stage is just because Ni somehow slows down the precipitation of 
Cu clusters. So Ni indeed has a temporal "delay" effect rather than refinement. 

It is then natural to apply the cluster mobility analysis for the ternary alloys. We simu- 
lated the mobility of Cu clusters in Fe-Cu-Ni alloys at 823K. As is shown in Fig. [HI diffusion 
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FIG. 7. Comparison of addition effect of third-party element on precipitation kinetics in view of 
mean size versus cluster density, (a) Ni effect in Fe-Cu-Ni alloys, (b) Zr effect in Al-Sc-Zr alloys. 

coefficient decreases linearly with area density of Ni, a minimum of about 50% lower dif- 
fusion coefficient than binary alloy is found for studied area densities. While the life time 
of clusters increases with area density and has a maximum of about 25% higher than bi- 
nary alloy for studied area densities. So the decreased diffusion coefficient is responsible for 
the delay effect of Ni. Even though, the slope of diffusion coefficient curves deceases with 
cluster size, larger clusters are less affected by the area density of Ni. So the delay effect 
by decreasing diffusion coefficient will be weaken with the clusters growing larger. Thus 
the Cu cluster number density of ternary alloys is expected to converge to a common value 
in the long term precipitation, which has been observed by experiments of Buswell et al^ 
and simulation of Bonny et air-. Also, the linear law of decreasing diffusion coefficient may 
explain the sequence of the delayed curves of advancement factor and number density of 
ternary alloys. We actually has an unexpected sequence of 1.5 at. % -0.5 at. % -1.0 at. 
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FIG. 8. Cluster mobility dependency of pA (Ni) at 823K. 

% -2.0 at. %, that is to say Fe-1.34 at. % Cu-1.5 at. % Ni is evolving faster than other 
ternary alloys in real time. While MC time gives an ordered 0.5 at. % -1.5 at. % -1.0 at. % 
-2.0 at. % sequence. Since has a hyperbolic relationship with Ni content, and diffusion 
coefficient has a linear relationship, numerically they will yield a minimum time on a specific 
Ni content, which means a ternary alloy containin this Ni content has fastest evolution. It 
is possible that this specific Ni content is near 1.5 at. %. 



V. CONCLUSION 

We have simulated the precipitation during thermal ageing of binary Fe-Cu and ternary 
Fe-Cu-Ni alloys by AKMC method. The energy model is based on a two-body short range 
model. A nonlinear time adjusting method from post-processing data has been proposed. 
Using the combined computational techniques, though the cluster density is over-estimated, 
good agreement of Cu precipitation kinetics has been obtained over a temperature range. 
For the effect of Ni on Cu precipitation, the following conclusions have been found: 

1. Peak number density of Cu clusters is higher in Fe-Cu-Ni ternary alloys, which can 
be explained by Ni promoting nucleation of Cu clusters. Surface area density of Ni on 
clusters has a linear relationship with Ni content; 

2. A delay effect has been found for Ni on Cu precipitation in the coarsening stage. 
Comparison with Al-Sc-Zr alloys reveals this effect is mainly temporal, rather than 
refinement of precipitations; 
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3. Diffusion coefficient deceases linearly with area density of Ni, while life time increases. 
This Deceasing effect of the cluster diffusion coefficient is responsible for the delay 
effect. Even though, this effect weakens with larger cluster size. 

On the other side, limitation of this AKMC study still remains. We are working on better 
activation energy model to over-come the shortcoming of pairwise energy model. 
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